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ABSTRACT 

Critically understanding the 'standard candlc'-likc behavior of Type la supernovae requires under- 
standing their explosion mechanism. One family of models for Type la Supernovae begins with a 
deflagration in a Carbon-Oxygen white dwarf which greatly accelerates through wrinkling and flame 
instabilities. While the planar speed and behavior of astrophysically-relevant flames is increasingly 
well understood, more complex behavior, such as the flame's response to stretch and curvature, has 
not been extensively explored in the astrophysical literature; this behavior can greatly enhance or 
suppress instabilities and local flame-wrinkling, which in turn can increase or decrease the bulk burn- 
ing rate. In this paper, we explore the effects of curvature on both nuclear flames and simpler model 
flames to understand the effect of curvature on the flame structure and speed. 

Subject headings: supernovae: general — white dwarfs - hydrodynamics — nuclear reactions, nucle- 
osynthesis, abundances — conduction — methods: numerical 
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1. INTRODUCTION 

Type la supernovae are used as 'standard candles' for 
cosm ology (see the review by iHillebrandt fc Niemeverl 
120001 and references therein). The standard model 
for a Type la supernova involves a flame that begins 
deep in the interior of a Carbon-Oxygen Chandrasckhar 
mass white dwarf, but the full mechanism for explo- 
sion i s not yet understood. One-dimensiona l simula- 
tions (jNomoto et al.lll984tlWooslev et al1l!984|) of flame- 
powered Type la can match the energetics of the obser- 
vations if the flame reaches ~l/3 of the speed of sound 
as it burns through the star, although the precise mecha- 
nism for this flame acceleration is still poorly understood. 
Other models begin with a flame that accelerates to 
-1/30 of the sound speed Ipommguez fc Hoflicbll2000|) 
and un dergoes an as-yet unexpla i ned transition to a deto- 
nation l|Arnett fc Livnell994albl : lKhokhlovll991albD . Ei- 
ther mode requires that the flame accelerate consider- 
ably beyond its laminar speed, and flame instabilities are 
generally cited as potential mechanisms. When burn- 
ing begins, pockets of buoyant, hot ash are formed as 
the flame propagates outward. This flame is unstable 
to the Landau-Darrieus and R ayleigh- Taylor instabilities 
ijBvchkov fc LibermanllT995aD . These instabilities wrin- 
kle the flame front, increasing the surface area, and there- 
fore, the bulk burning rate. 

Both strain from curvature and flow-induced stretch 
of a f lame are known t o affect a flame's speed and struc- 
ture (|Marksteinl Il964|) . These strains change both the 
local burning rate, and the bulk burning behavior in com- 
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plex large scale flows (|Helenbrook fc Lawlll998l) through 
multidimensional flame instabilities. Curvature effects in 
terrestrial flames have received a gre at amount of atten - 
tion (see, for instance, the review bv lLaw fc Sunell2000fl . 
but the astrophysical flame literature on these effects has 
been sparse. 

The instabilities that might increase the total burn- 
ing by a flame do so by stretching and curving the flame. 
Because these strains themselves can significantly change 
a flame's local burning rate, the onset and growth rate 
of flame instabilities are directly affected by curvature. 
Thus, in the context of the standard Type la model, un- 
derstanding the detailed explosion mechanism requires 
understanding the micro-physics of the flame propaga- 
tion in the white dwarf. One way of doin g this would be 
to ex tend large supern ova si mulations feg~ |Gamezo el; alJ 
120031 INiemever et allll99fit IReinecke et al.lll999L 12002b) 
to incorporate a fully resolved flame. However, at the 
onset of the burning, the flame is only ~ 10 -5 cm thick, 
inside a star of radius ~ 10 s cm. At lower densities, the 
flame becomes broader, but remains tiny compared to the 
radius of the whole star. Thus any realistic study of mul- 
tidimensional Type la scenarios that follow a significant 
portion of the star requires understanding of the behavior 
of the flames, and must use a model on unresolved scales 
to describe the burning physics, such as flame-stretch in- 
teractions and the corre sponding effects on instabilities- 
Some such models exist l)Khokhlovlll995t IReinecke et alJ 
I2002a|t . but they are based on the properties of large 
scale flows and terrestrial scalings, rather than on the 
local physical properties of astrophysical flames. A de- 
tailed understanding of the flame's response to curvature 
is necessary input to a subgrid scale flame model. 

In this paper, we describe numerical experiments 
of propagating, spherical flames, similar to re- 
cent experiments and calculations performed b y the 
chemical combustion community iAmi^e^^^ 1997: 
I Bradlev et all 1199(1 Il998t IGu et alJ 1200(1 IHassan et alJ 
119981 IKarnov et al.111 99ft IT" 



H998t " lKarDOv ^rSnil997f ISiin et al.l l1999tl. designed to 
measure the flame structure and speed response to geo- 
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metric curvature. All of the flame simulations presented 
are 1-d (either in planar or spherical geometry), to get 
clean and unambiguous measurements of the flame's re- 
sponse to curvature. As the flame propagates to different 
radii, it experiences different strain due to the geometry, 
parameterized by a dimcnsionlcss strain rate, Ka, defined 
in later sections. Quantitatively, we seek the modified 
flame velocity 

S t = S° ■ /(Ka), (1) 

where Si is the laminar flame speed, and Sf is the un- 
curved laminar flame speed. (Throughout this paper, 
a superscript will refer to planar, unstretched flame 
q uantities.) 

iMarksteinl l|1964D postulated a flame-speed law of the 
above form for small curvatures based on empirical re- 
sults. Expressed in dimcnsionlcss numbers, this is writ- 
ten: 

S { = 5j°-(l + MaKa), (2) 

where Ma is the dimensionless Markstein number, which 
depends on parameters of the flame. We refer to this 
as the Markstein relation. The Markstein number can 
be negative or positive for any given flame. A negative 
Markstein number indicates that the flame would burn 
more slowly in regions of positive curvature (Ka > 0) and 
more quickly in regions of negative curvature (Ka < 0), 
and the reverse is true for a flame with positive Markstein 
number. 

The Markstein relation has proven very robust for de- 
scribing flame behavior in terrestrial flames. In that con- 
text, work has been done to calculate Markstein num- 
bers f rom first principles (for example , iChung fc Lawl 
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In this paper, we discuss its applicability to astrophysical 
thermonuclear flames and examine the flame-curvature 
relation both for these flames and for simpler model 
flames with parameters that may be varied to understand 
their effects. 

We outline some theory relevant to our study in [J2 
In ^Jwe discuss the methods used in the numerical ex- 
periments. These experiments, which are described in 
21 were performed to understand how the flame speed 
changes with curvature. In fJSlwc present the results of 
these experiments, and we conclude in !£2 The appendix 
provides some simple tests that provide verification of 
our code for the present simulations. 

2. THEORY 

2.1. Planar Astrophysical Flame Structure 

The structure of an astrophysical flame is sketched in 
Fig. ^ Between the unburncd state (which through- 
out this paper we will label with a subscript u) and the 
burned state (subscript 6), there are two distinct regions. 
There is a 'preheat' zone, where energy from the reac- 
tions and the already-hot ash diffuses outwards and heats 
the incoming (in the frame of the flame) fuel, and there is 
a reaction zone where the bulk of the thermonuclear reac- 
tions occur. The interface between these two zones — the 
state at which burning begins — we will denote with sub- 
script i. Because the flame moves very slowly compared 
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Fig. 1. — A schematic of an astrophysical thermonuclear flame, 
propagating to the right. The exothermic reactions arc largely 
concentrated in a very thin reaction zone, and energy from that 
region diffuses into the cold fuel over a large preheat zone. 



to th e sound speed in the domain ijTimmes fc Wooslevl 
Il992|) . the pressure is nearly constant throughout the un- 
burned fuel and burned ash, Pb ~ P u = P ■ It is possible, 
however, for the pressure to differ significantly from this 
value in the flame structure itself, because of the exother- 
mic reactions taking place there. 

The astrophysical flames of interest to us here are 'pre- 
mixed' — that is, the combustible fuel is waiting only for 
sufficient heat to ignite. Thus, the speed of propagation 
of the flame is determined completely by two rates - 
the rate of energy input by the reaction, and the rate 
at which enough heat is diffused outwards to ignite the 
as- yet unburncd fuel. 

Astrophysical premixed flames differ from their ter- 
restrial counterparts in several respects. If we were to 
sketch a terrestrial flame in Fig. ^ we would have, over- 
lapping with the thermal diffusion zone, a material dif- 
fusion zone where fuel diffuses inwards towards the fuel- 
depleted burned zone. The Lewis number, Le, describes 
the relative importance of thermal and species transport 
across a flame; in terms of the thermal diffusivity D t h 
and the material diffusivity D Sl Le = D t h/D s . In the 
astrophysical flames we consider here, the Le wis num- 
ber is of order 10 7 <)Timmes fc Wooslevlll992h . whereas 
in terrestrial flames, the Lewis number is often of or- 
der unity. Since for astrophysical flames, the effect of 
species diffusion is insignificant compared with thermal 
diffusion, we neglect species diffusion for astrophysical 
and some model flames; in those cases, we informally 
speak of an infinite Lewis number. Note, however, that 
in these numerical simulations, numerical diffusion will 
result in some small species diffusion. 

Astrophysical flames also have extremely peaked 
energy-generation rates - - that is, rates which are 
strongly dependent on temperature. This is usually de- 
scribed by a Zeldovich number, Ze, defined as 



Zc 



T a Tb — T u 



(3) 



Here, T a is an activation temperature that sets the burn- 
ing rate by representing an activation energy for the re- 
action, such as a potential energy barrier that must be 
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overcome for the reaction to proceed, and T\, and T u are 
the temperatures of the burned and unburned gas. In 
hydroc arbon-air burnin g, this number is typically of or- 
der 10 JgI assmanlfl996|) : most of the flames we will see 
here will have similar ratios, meaning similarly peaked 
energy generation rates. 

Astrophysical flames also typically occur in electron- 
degenerate material, meaning that the large amounts of 
energy input to the flow (much larger compared to the 
ambient energy density than in terrestrial flames) result 
in comparatively small density changes at the roughly 
constant pressure relevant for flame problems. The ra- 
tio of the burned-gas density to the unburned-gas density 
will be parameterized by a = p u /pt- Another factor con- 
tributing to a relatively low a in astrophysical flames is 
that these flames arc powered by fusion, so that the mean 
molecular weight of the ash behind the flame is, unlike 
many chemical flames, greater than the mean molecular 
weight of the fuel ahead of it . 

The reactions arc governed by a strongly-peaked burn- 
ing rate (ie., large Ze), so that in the preheat zone, the 
fuel can get quite hot without any significant burning 
taking place; instead all of the burning occurs only in the 
hottest region, in a narrow burning zone. In this case, 
a metho d for estimating the flame veloci t y whic h dates 
back to iZeldovich fc Frank-Kamenetskiil l|1938j) shows 
the most important physical effects. For a planar, steady 
flame, one can consider the energy equation in the refer- 
ence frame of the flame front itself, and match boundary 
conditions between the non-burning preheat zone where 
thermal conductivity is important, and a narrow reac- 
tion zone. If one uses a canonical peaked temperature 
dependence for a burning rate used in combustion the- 
or y, an Arrhenius la w (see, for e xampl e, textbooks such 
as iGlassmanl l)1996j) : iWilliamsl l)1985ft : IZeldovich et~aTI 
l| 19851) ). which is exponential in the ratio —T a /T, and 
consider our nuclear generation rate to be from a second- 
order (eg., fuel + fuel — » products, so that the reaction 
depends on the square of the fuel concentration) then the 
reaction rate is 



apXje 



2 P ~TJT 



(4) 



where Xt is the fuel mass fraction, then one can find 



S? 



Oh 



\c P f b 



-T a /T b 



(5) 



where a is the thermal conductivity and C p is the specific 
heat at constant pressure. 

This approach is sufficient to get velocities for the as- 
trophysical flames we describe here to within a factor of 
a few, and demonstrates the most important physical ef- 
fects. A more careful analysis bv IBychkov fc LibermanI 
l)1995b|) . using asymptotics and explicitly including an 
equation of state and conductivity suitable for com- 
pletely degenerate materials, can model the laminar 
flame speeds to within 25-50%. 

2.2. Wrinkling 

Above, we have assumed a planar flame. Now con- 
sider Fig. |3 In curved regions of the flame, heat trans- 
port into the cold fuel can be concentrated or diluted 
by simple geometry of the wrinkling, affecting the flame 
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Fig. 2. — A schematic of a wrinkled flame, with the direction 
of propagation shown by the large arrow. Heat diffuses from the 
burning zone to the cool fuel as shown by the little arrows across the 
flame. At regions of positive (negative) curvature, heat diffusion 
is 'diluted' (concentrated), and flame propagation will slow down 
(speed up). 



speed and structure in those regions. Since astrophys- 
ical flames propagate by thermal conduction, when the 
heat transport is concentrated by negative curvature, we 
would expect the flame to propagate faster, and when 
heat transport is diluted by positive curvature, we would 
expect the flame to propagate more slowly. 

Note that this figure does not describe cusping of an 
interface moving at constant velocity, such as in the 
Landau-Darrieus instability. Here, we emphasize that 
the flame speed itself will change along the flame; if the 
flame speed were to vary strongly enough, cusping need 
not occur. 

The change in geometry can be taken into account by 
including area terms into the equations of the previous 
section, and we can derive results that include geometric 
eff ects due to curvatu re. This was done, for instance, 
in lLaw fc Sund lj2000|) . Defining the strain due to the 
curvature with a dimcnsionlcss Karlovitz number, 



Ka c 



dA 
~A 



A(x u ) - A{x b ) 
A{x u ) 



1 dA 
A~dx' 



(0) 



where Xb and x u are the position of the burned and un- 
burned material, respectively, and If is the thickness of 
the flame. We then find to linear order in Kap, 



S? 1 



and so for our large-Lewis number flames, 



-^Ka c 



(7) 



(8) 



Note that the Karlovitz number is defined so that it is 
negative if the area at the burned fluid is greater than 
that at the unburned fluid, ie. for regions of negative 
curvature. Thus, since Ze > 0, we get flame speed en- 
hancement in negatively curved regions, which is as we 
expect. Also notice that the effect is strongest for flames 
with very peaked temperature dependencies, where the 
geometrically-induced difference in temperature diffusion 
will have a large effect. 
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2.3. Stretch due to flow 

A steady spherical flame will feel the strain described 
above. Another contribution to strain in the flame comes 
from the flow — for instance, a flame propagating out- 
ward will also be advected outward by the expansion of 
material behind it. Any such flow will also lead to a 
strain, characterized by a strain rate K = (l/A)(dA/dt) 
where A is the area of the flame surface. For the spherical 
flames we consider here, A = Airr 2 , so that 



K = — u 



gas* 



(9) 



where rf is the radius of the flame's position, and it gas 
is the flow velocity of the gas at the flame's location. 
This strain can be converted to a dimcnsionless Karlovitz 
number by considering the strain over a flame crossing 
time, lf/S$: 

Ka s = 2 ±^. (10) 
r S? 

Note that for a flame in a medium where the flow is zero, 
this strain vanishes. 

One can also consider a total strain rate, Kay = 
Kac+Kag. We can meaningfully do this for the purposes 
of ass essing the effect on flame velocity; IGroot et al.l 
l)2002|) demonstrated that a single Markstein number can 
be found that is meaningful for the combined effects of 
stretch and curvature, and in fact considering only the 
stretch is problematic. We note here that other works 
in the literature choose different sign conventions for Ka 
and Ma. 

2.4. Effect on Landau-Darrieus growth rate 

The change in flame speed as a function of curvature 
can modify multidimension al flame instabilities. T he 
Landau-Darrieus instability ijLandau fc Lifshit3ll987tl is 
often proposed as a mechanism for accelerating the burn- 
ing rate of a flame in a white dwarf. In an otherwise 
planar flame with a density ratio a = p u j pb, a sinusoidal 
'wrinkle' of wavenumber k will grow in the linear regime 
with an exponential growth rate 



(11) 




a + 1 



One can re-derive the linear theory of Landau-Darrieus 
growth rates with the non -constant flame speed given in 
Equation [3 and arrive at l)Zeldovich et al.lll985j) 



burn more slowly, because they're at positive curvature, 
so they tend to fall back, and the 'troughs' of the flame 
at negative curvature burn faster and catch back up. 

For linear perturbations, the increased burning rate 
caused by a linear curvature term at 'troughs' in the 
wrinkle is exactly offset by a decreased burning rate at 
the 'peaks', so that the total burning rate remains pro- 
portional to the surface area of the flame, as in the con- 
stant flame velocity case. Thus the total burning rate 
is decreased by the decreased growth rate of the insta- 
bility. If cusping occurs when the perturbation becomes 
nonlinear, then (fast-burning) regions of negative curva- 
ture decrease to almost zero in favor of regions of (slow- 
burning) positive curvature, which further decreases the 
net burning rate. 

Thus, for understanding the behavior of the Landau- 
Darrieus instability in astrophysical phenomenon, at 
least the curvature behavior of the flame must be mod- 
eled. An alternative is to model the ent ire flame behav- 
ior, as in IBvchkov fc LibermanI l)1995al) . where a care- 
ful analysis of the linear growth rate and linear in- 
stability of a degenerate thermonuclear flame to the 
Landau-Darrieus instability was performed. These au- 
thors also find a cut- off beyond which modes a re sta - 
ble. Another work, INiemever fc Hillebrandti l)1995j) . 
studied the instability numerically, although that work 
obtained anomalous results. In both analytic cases, 
there is a most-uns table mode; there is some evidence 

1mm 



i|Helenbrook fc Law! 119981) that it is this mode which 
grows fastest in a flame going through a complex or tur- 
bulent flow field if there is p ower on this scale, and indeed 
there are other suggestions l)Denetill998j) that small scale 
wrinkling can directly induce larger-scale wrinkling. 

3. NUMERICAL METHODS 

We consider two classes of flames — 'model' flames and 
'astrophysical' flames. The delineation is reflected in the 
input physics. In the simpler, 'model' case, we use a poly- 
nomial or simple exponential expression for the burning 
and a simple ideal gas equation of state. At the complex, 
astrophysical end, we use a relativistic, degenerate equa- 
tion of state (EOS) and a more astrophysically relevant 
reaction rate. In this section we review the physics used 
in the astrophysical-flame simulations. Similar details for 
the model flames are described in Appendix |0 




kl f M&(kl /Ma + 2a) - 1 + klj 
(12) 

where If is the therm al width of the flame. In fact, some 
authors (for example . iClanet fc Searbvl (11998^ 1 have used 
this to measure the Markstein number of a flame from an 
experimentally determined growth rate of the instability. 

From our discussion of the effects of positive and neg- 
ative Markstein number, we expect our astrophysical 
flames to have Ma < 0. In this case, the Landau- 
Darrieus growth rate is reduced, and indeed modes 
k > (a — l)/(2a^|Ma|) are stable against the Landau- 
Darrieus instability. In a negative Markstein number 
flame, the flame-curvature relation works against a wrin- 
kle trying to increase in size; the 'peaks' of the flame 



3.1. Hydrodynamics 

For the hydrody n amics , we use the PPM 
ijColella fc WoodwardI Il984|) co mpressible hy drody- 
namics module in the Flash code (|Frvxellll200(i(i . PPM 
is/ a widely used Godunov method, which solves the Eu- 
ler equations in conservative form using a finite-volume 
discretizatio n. Our impleme ntation of PPM is described 
in depth in iFrvxelll l|2000D . and the P PM module i n 
particular has been tested rigorously in iCalderl l)2002j) . 
The Riemann solver is capable of dealing with a general 
equation of state, f ollowing the method outlined in 
IColellafc Gialfmsl . 

To accurately follow reactive flows, a separate advec- 
tion equation is solved for each species, 

^i + ^.pX i v = pR i , (13) 

where Xj is the mass fraction of species i, subject to the 
constraint that they sum to unity, and Ri is the net rate 
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of change in the abundance of species i due to the nuclear 
burning. e nuc and Ri are related by 

£nuc = 51 e&i-Ri , (14) 
t 

where tu is the binding energy per mass of species i. 

3.2. Equation Of State 

Astrophysical flames in degenerate white dwarf mate- 
rial have a complicated equation of state. Contributions 
from the ions, electrons, and radiation must be accounted 
for. The ions are assumed to be fully ionized, and we 
use an ideal gas expression for them, as described above. 
Finally, the radiation term is simply blackbody. Full dc- 
tails of the implementati on of this EOS can be found in 
iTimmes fc Swestvl 

The equation of state is responsible for much of the 
character of astrophysical flames. For the conditions that 
we consider, the pressure is dominated by the degenerate 
electron contribution. The degeneracy means that pres- 
sure responds only weakly to temperature changes, so 
the resulting density jumps behind the flame are smaller 
than they would be with a pure ideal gas EOS (and ev- 
erything else kept the same). 

3.3. Diffusion 

Diffusion is implemented in Flash with explicit time 
differencing in an operator-split manner. The heat flux 
term is added to the energy evolution equation 

dpE 



dt 



V • (pE + P) v = V • (crVT) + (Xi,p,T). 



(15) 

Here, E is the total energy per unit mass, P is the pres- 
sure, v is the velocity, T is the temperature, <r is the 
conductivity, e nuc is the energy source term from burn- 
ing, X t is the mass fraction of species i, and p is the 
total density. The heat flux — cVT is computed during 
the hydrodynamic step and added to the total energy flux 
computed by the PPM solver. The fluxes are then used 
to update the energy in each zone. This implementation 
preserves the conservative nature of the algorithm. 

Because the diffusion is calculated explicitly, the diffu- 
sion term adds an additional timestep constraint 



1 Sx 2 

dt < 

" 2D th 



(16) 



where D t h is the diffusion coefficient, D t h = a/(pCp). 
The simulation evolves at either the diffusion timestep 
or the CFL timestep — whichever is smaller. 

For our more realistic thermonuclear flames, the ther- 
mal diffusion is the only diffusive process modelled. For 
some of our model flames, we allow finite Le, and so 
species diffusion is also calculated; the computation pro- 
ceeds in the same manner. This is described in Ap- 
pendix [nj 

Since these diffusion mod ules were not included in the 
original Flash code paper l)Frvxel]ll20f)Ot ). we include in 
Appendix ^ tests of the correctness of the diffusion op- 
erators in Cartesian and spherical coordinates. 

3.4. Conductivity 

For testing purposes and for our model flame propaga- 
tion problems in Appendix IO we use a simple constant 



diffusivity to describe the heat transfer. However, the 
conductivity for astrophysical flames of interest to the 
Type la problem is much more complicated. For degen- 
erate carbon, electron-electron and electron-ion collisions 
are important processes in heat transport. The conduc- 
tivity in the post-flame state can be some 3 orders of 
magnitude higher than that of the pre-flame state, mak- 
ing a constant conductivity a bad approximation. We 
use a stellar conductivity routine th at includes thes e pro- 
cesses as well as radiative opacities l)Timmesll200Cl() . 

3.5. Burning 

We experiment with several reaction rates and net- 
works for our flame simulations. In this paper, we use 
one-step irreversible reactions, with increasingly complex 
reaction rates. In all cases the energy release from the 
nuclear burning is computed in an operator split fash- 
ion. After the hydrodynamics are evolved, the reaction 
network computes the energy release and change in nu- 
clear abundances over the course of the hydro timestep. 
This energy release is then added to the internal en- 
ergy of that zone, and a new temperature is computed. 
For these fully-resolved flames, the CFL-limitcd timestep 
from the explicit hydrodynamics is much smaller than 
the timescale to significantly change the temperature by 
burning, so that explicit coupling is not needed. We de- 
scribe the burning rates for our astrophysical flames here; 
model flame burning is described in Appendix IU1 

3.5.1. One-step Carbon Burning 

The astrophysically-motivated reaction we use is a 
one-step 12 C + 12 C reac tion. The rate comes from 
iCaughlan fc Fowlerl l|1988ft (we refer to this as the CF88 
reaction network) and converts two 12 C nuclei into a sin- 
gle 24 Mg nuclei , releasing 5.57 x lQ - ^erg g " 1 . Following 
the notation of ICaughlan fc Fowlerl l)1988[) . the reaction 
rate is 

" = ~XhpN A X, (17) 



dt 



where 



,T, 



5/6 



N A X = 4.27xI0 26 ^ T exp \ 84165 - 2.12 x lO' 3 ^ 



T, 



3/2 



T, 



1/3 



(l,o 



(18) 

where T 9 = T/10 9 K, T 9iQ = T 9 /1.0 + 0.0396T 9 , and 
Xc is the mass fraction of 12 C. For each 24 Mg nucleus 
created, 13.933 MeV is released. 

This is the main reaction involved in a pure carbon 
flame in a white dwarf. We treat all other species (eg. 
16 0) simply as inert dilutants that serve only to reduce 
the effective burning rate. We assume here that other re- 
actions (eg. 24 Mg(a,7) 28 Si burning), if modeled, would 
occur well behind the flame and would not be as impor- 
tant in setting the properties of the flame. 

With this reaction mechanism, we neglect any ion 
screening, for simplicity. Screening is not an impor- 
tant effect except at low temperatures or high densities. 
At high densities, a one-stage reaction as modeled by 
this rate might have other problems; late stage reactions 
could be important. Conversely, low temperatures are 
not of interest to us for these flame-propagation prob- 
lems. Thus we restrict ourselves to regimes where screen- 
ing is unimportant. 
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3.6. Boundary Conditions 

Boundary conditions for subsonic flows in a compress- 
ible code arc non-trivial. The zero-gradient /outflow con- 
dition usually used is insufficient, since not all the charac- 
teristics at the edge of the computational domain will be 
leaving the domain. This means that some information 
can enter the domain from the boundary and introduce 
noise into the flame solution. 

Boundary conditions for a subsonic outflow have been 
propo sed (see for example iHedstroml Il979t iThompsonl 
11987(1 . but no well-accepted general solutions exist. 
These boundary conditions try to extend the solution 
into the boundary conditions in such a way that informa- 
tion is not transmitted into the domain from the bound- 
aries. We choose problem-specific boundary conditions 
that allow most waves to leave the domain. We first 
use the fact that, to a good approximation, the pres- 
sure is constant in the burned and unburned regions, 
and since the flames are not in a closed vessel, the pres- 
sures will be constant in time. This allows us to set the 
pressure in the boundary conditions as a constant. Fur- 
thermore, since we want heat to diffuse out, the tempera- 
ture boundary condition is simply set to be zero-gradient. 
The boundary conditions for the mass fractions are also 
zero-gradient. This gives us enough information to find a 
thcrmodynamically consistent density through our equa- 
tion of state. Finally, the velocity in the boundary con- 
ditions is found by enforcing mass conservation, 

p l v i A i = p i+ iv i+ iA i+ i (19) 

if the velocities arc out of the domain, or zero if the veloc- 
ities are inwards. Here Ai is the area of the face through 
which the material flows; in planar geometry, this cancels 
out, but it is important in the spherical geometries we 
consider. A demonstration that this boundary condition 
functions in the way we'd like is given in Appendix 1X1 

For our spherical runs, we use at r = a reflect- 
ing boundary condition. This boundary condition is re- 
quired by the geometry; spherical symmetry requires, for 
any quantity f(r) other than the radial velocity, f(r) = 
/(— r), and for the radial velocity, v r (r) — —v r {—r). This 
is exactly the reflecting bo undary condi tion as imple- 
mented in the Flash code l|Frvxellll20"00(l . 

3.7. Adaptive Mesh Refinement 

Flash employs an adaptive mesh thr ough the 
PARAMESH package (jMacNeice et al J 12000(1 . allowing 
it to put computational zones where the resolution is 
needed and to follow smooth flow with less resolution. 
The Flash mesh adapts by dividing the domain in half 
in each dimension, creating new blocks (2 in 1-d, 4 in 
2-d, 8 in 3-d) that are logically the same as their parent 
(typically containing 8 computational zones in each di- 
rection) but with twice the spatial resolution. Each new 
block is checked to see if more refinement is needed, and 
if so, it is further subdivided. When a block is refined, 
the newly created children need to be initialized with 
information from their parent. To obtain higher-order 
accuracy, a quadratic polynomial is fit to the data of the 
parent zone and the two zones on either side of it. This 
procedure is described in detail in Appendix iBl 

When evolving a flame, it is only necessary to put the 
resolution near the reaction and diffusion zones. We re- 
fine when the second derivatives of density or pressure 



and are large (the former is expected to be large at the 
fuel/ash interface, and the latter in the burning region). 
We also refine on the second derivatives of the nuclear 
energy generation rate and temperature. The methods 
for evaluating the derivatives a nd determinin g their mag- 
nitude is discussed elsewhere l|Frvxellll2000(l . These cri- 
teria ensure that the region around the flame is highly 
resolved. 

4. EXPERIMENTS 

4.1. Overview 

To numerically investigate the effect of curvature on 
astrophysical flame speed, we conduct numerical ex- 
periments of 1-d model and astrophysical flames run- 
ning out from, and in toward, the center of a spher- 
ical domain, as has been done in the terrestrial com- 
bustion c ommunity experimentally and compu t ation- 
ally (eg. lAune et all Il99& iBradlev et alJ Il995 pSl 
Gu et all 120001 lHassan et al.lll998t iKarpov et all 119971 
Kwon fe FaetU200lUMuller et al.ll997tlSun et al.ll999() . 
Doing so allows us to examine the flame speed and struc- 
ture for a range of curvatures and strain rates, both pos- 
itive and negative, by taking snapshots of the flame at 
varying positions (and thus curvatures) throughout the 
domain. We simulate the same flames in planar geometry 
for comparison. 

Focusing on 1-d is computationally more efficient and 
allows us to separate the pure strain and curvature ef- 
fects from more complicated multidimensional instabili- 
ties, which will be the focus of a later paper. In addition, 
for the present work, we studied the effects of resolution 
on the experiments a nd m easurements. The details arc 
discussed in Appendix lA.4l In the numerical experiments 
we present, the flames were sufficiently well-resolved to 
accurately model the flame structure and speed. 

We begin with the simpler model flames, adjusting pa- 
rameters to ensure we understand the effects of each. For 
each set of input physics (EOS, reaction mechanism, and 
model parameters), we simulate three flames — a planar 
flame, a spherical flame propagating out from the center, 
and a spherical flame propagating toward the center. We 
then continue on with simulations of a simplified astro- 
physical reaction rate with a complex EOS. 

In each case, we ignite the flame by placing hot ash 
- at a temperature corresponding to the ambient pla- 
nar burned temperature of the flame, unless stated oth- 
erwise — alongside the cold fuel we wish to burn, in 
pressure equilibrium. Neglecting to place the fuel and 
ash in pressure equilibrium will generate strong pressure 
waves, especially in non- or partially degenerate fluids, 
giving spurious results — one is no longer studying a 
flame, but a transient pressure-driven reaction front. We 
be lieve this is the reason for the a nomalous results found 
in lNiemever fc HillebrandU (|1995tl . which have not been 
reproduced in work since. 

For planar flames, and for the outer boundary con- 
ditions for our flames in spherical coordinates, we use 
the constant-pressure boundary conditions described in 
For the inner boundary conditions in spherical co- 
ordinates, we use reflecting (or 'symmetry') boundary 
conditions at x = to be consistent with the spheri- 
cal geometry. These boundary conditions are consistent 
with those used in the terrestrial combustion literature, 
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Fig. 3. — Position versus time, measured four different ways, 
for both a KPP (left) and an Arrhenius (right) flame. Both flames 
have Le = oo and are propagating spherically outwards. Results 
from all four position measures are plotted on the top. On the 
bottom arc plotted differences from using {e n ucx) / (e n uc). The 
three other methods are: maximum e n uc (+), maximum d 2 T/dx 2 , 
(-), and maximum dXf/dx (o). The KPP flame has a reaction zone 
of 6 cm, and the Arrhenius flame has a reaction zone of 1/2 cm. 
The Arrhenius reaction takes about 50 seconds to ignite a self- 
sustaining flame; before then, any measure of 'flame position' is 
necessarily meaningless. 



eg., l)Sun et al.lli999j) . 

4.2. Measuring flame position 

One of the most fundamental things we need to do is 
to measure the flame's position (r/) at different times. 
Since flames have finit e thickness, 'the' fla me position is 
ambiguous. Following lOroot ~ we take the 

name's position to be that of the reaction zone, which 
is less sensitive to strain and curvature effects than, 
for instance, the outside of the preheat zone. Then a 
consistent Markstein number can be calculated that in- 
cludes the effects of curvature and stretch. Since even 
the reaction zone has finite thickness, we consider sev- 
eral measures for the position of the reaction zone. The 
first is the most obvious, the position of maximum e nuc . 
This has the advantage of simplicity, but could poten- 
tially be noisy, since extrema are generally not numeri- 
cally well behaved. Another is (e nuc a;) / (e nuc ), an e nuc - 
weighted position. Since e nuc is very strongly peaked, 
this should localize the flame. Clearly, in the limit of a 
delta- function e nuc , this measure reduces to finding the 
x of maximum e nuc , but as an integral quantity, it will 
vary more smoothly. A third measure is where the fuel 
concentration is most rapidly changing — the location of 
maximum dXf/dx. A fourth is a temperature criterion, 
the location of maximum d 2 T/dx 2 . In Fig. we show 
the results of using these measures for a simulation with 
a KPP flame and an Arrhenius flame. The KPP flame 
has a reaction zone approximately 6 cm wide in a domain 
of 500 cm, and the Arrhenius flame has a reaction zone 
of approximately 1/2 cm in a domain of 100 cm. 

Because the KPP reaction zone is so wide and includes 
the preheat zone, finding an unambiguous 'location' is es- 
pecially difficult. The reaction rate depends as strongly 
on fuel as on temperature, so the fuel-change rate is 
fastest right at the front of the flame. Also, the energy 
generation rate is large in a wide region, so averaging 
over e nuc produces a position that lags well behind the 



front of the zone. 

For the Arrhenius rate, however, once ignition occurs, 
all of the direct burning measures (e n uc, (a^nuc) / (enuc); 
and dXf/dx) fall almost exactly on top of each other, 
with only the temperature measure leading slightly, 
as the temperature must start diffusing outward from 
the burning region to change the temperature's second 
derivative. 

In light of the above, we use the e nuc - weighted posi- 
tion as a marker for the position of the reaction zone, 
and thus the flame, for flames with strongly-peaked re- 
action networks (Arrhenius, CF88) as it generates fairly 
clean velocities through differencing and is otherwise es- 
sentially degenerate with the very intuitive maximum-of- 
£nuc measure of position. For KPP, we use the position 
of maximum e nuc . 

4.3. Measuring flame thickness 

The thermal thickness of the flame is an important 
quantity, setting the length scale for relevant microphys- 
ical effects. Since this width is set by a diffusive process, 
however, its precise start and end arc poorly defined. We 
use two measures of the flame thickness, both used in the 
literature, that span the range of reasonable measures of 
the thickness. Both are based on the change in tempera- 
ture from the burned to unburned state, AT = Tb — T u . 

The first way we use to measure flame thickness, 
method I, is the '10/90' a pproach used for instance in 
iTimmes Wooslevl lTl99^ . It measures the distance 
from where the temperature crosses 10% above the un- 
burned temperature to 90% of the way to the burned 
temperature 

lf> = x(T = T b - .1AT) - x(T = T U + .1AT). (20) 

This gives quite wide thermal widths, as it measures well 
into the diffusive tails of the thermal structure. 

The second m easure, method II, used for instance in 
iSun et aP 1)1999(1 . finds the width by dividing the temper- 
ature difference by the maximum temperature gradient 
in the flame 

J max | Vi | 

This measure gives the thinnest meaningful thermal 
thickness of the flame, as any smaller thickness must im- 
ply a larger thermal gradient than any that exists in the 
flame structure. We report both flame thicknesses and 
use them to constrain measurement uncertainties in de- 
rived quantities that depend on flame thickness. 

4.4. Measuring flame velocity 

The velocity we are interested in is not the change of 
flame position over time (drf/dt), but the rate at which 
the flame is consuming fuel. This means we want to 
measure the flame speed with respect to the fuel. The 
expansion of hot burned ash smoothly generates a ve- 
locity field as the flame passes through. In the case of 
spherical geometry, the velocity of the material between 
the flame and the origin is zero, since the material has 
no place to go, and thus the only nonzero velocities are 
exterior to the flame. In the open planar geometry, ma- 
terial can flow freely through both boundaries. In both 
cases, the continuity condition requires Apu = const on 
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Fig. 4. — The gas velocity structure of a propagating flame 
with the flame propagating outward through a spherical domain 
(top), through a planar domain (middle), and inward in a spher- 
ical domain (bottom). In the closed geometries, the left state, 
cither burned or unburned, has a velocity of zero imposed by the 
boundary and the continuity condition, Apv = const. 



either side of the flame. Fig. ^sketches idealized velocity 
profiles for the three configurations. 

To find the propagation velocity into the fuel, one ap- 
proach is to find the change in the flame's position over 
time (by differencing, or by fitting and then differentiat- 
ing) and to subtract off the gas velocity (which can be 
directly measured). This can be problematic, however, 
especially in the outward-expanding flame case. The ex- 
panding gas, especially for large density contrasts, will be 
a large component of drf/dt, and thus any subtraction 
is likely to be noisy — especially since the flame velocity 
and the gas velocity are found in different ways. The 
situation is worse in the spherical case, as gas velocity 
has a non-constant spatial structure, adding difficulty to 
measuring 'the' gas velocity. 

Another approach, used in the chem ical literature (eg., 
lAung et aDll997t lHassan et al.|ll998j) . is to note that if 
the flame is a discontinuity, it must separate two regions 
with the same momentum flux in the frame of the flame. 
Given that in the lab frame, the velocity interior to the 
flame is zero, the flame propagation velocity must be 
{Pb/pu){drf /dt). However, this is an approximation that 
assumes an infinitely thin flame, and we are measuring ef- 
fects explicitly due to the finite structure of the flame. A 
better approximation, {A^pb) / {A u p u )[dr f / dt), includes 
the effects of geometry but requires locating 'the' posi- 
tion of the unburned and burned states, adding uncer- 
tainty due to the ambiguity of those locations. 

Since what we are really interested in is the burning 
rate, which is directly measurable, another approach is 
simply to measure the total burning rate, E = J e nuc pdV . 
If one then assumes that the reaction region is thin (a 
better approximation than assuming the thermal struc- 
ture of the flame to be thin), then one can compute 
a per-area burning rate, E / '(Airr'j). Knowing the per- 
area burning rate and the chemical energy per volume 




Fig. 5. — The propagation velocity - the speed of the flame 
into the fuel - of a Le = 5, Ze « 5 Arrhcnius flame propagating 

outwards in a spherical domain. Shown is (dr f /dt)(p b / p u ) (' '), 

{A b p b )/(A n p u )(drj /dt) (' — '), velocity from total e nuc ('o'), and 
(dr f /dt) - v gas ('+'). 



(AepXf/[if) of the unburned state, where Ae is the bind- 
ing energy per fuel particle and /i/ is its mass, one can 
find the speed at which the flame must be propagating 
through the unburned state to generate the measured 
energy, E p / / (Airr'j pX / Ae) . 

A comparison of several methods of measuring ve- 
locity is given in Fig. [SJ These plots were made 
for a typical flame with typically noisy velocity data, 
(the Le = 5, Ze = 5 outwardly propagating Arrhe- 
nius flame) showing the flame's propagation velocity 
through the fuel as a function of time; ignition occurs 
at t rs 50 sec. To calculate the area-corrected ex- 
pression (Abpt) j (A u p u )(dr f I dt) , we (rather crudely) as- 
sumed A b pa 4ir{r f -l° f {I) /2) 2 and A u « A-K{r f + f f (I) /2) 2 . 
Note that ignoring the thickness of the flame spuri- 
ously increases the measured effect of curvature in the 
(pb/Pu) {drf/dt) velocity measure. 

Another thing to note in Fig. [5] is that in the di- 
rect velocity measure drf/dt, one can see pulsations 
due to an essentially translational 1-d instability in 
a flame with large Lewis n u mber (see, for i nstanc e, 
iBvchkov fc Libermanl l)1995b|) : iZeldovich et~aTI Ip^ . 
While this slightly affects the instantaneous position of 
the flame front, and, therefore, velocity measures based 
on differencing the flame position, it has little effect on 
the total burning occurring in the flame, and thus in ve- 
locity measures based on bulk burning rate. 

4.5. Igniting the flame 

Our flames are ignited by placing a hot region of ash 
next to the fuel, and letting the heat diffuse into the fuel 
until the fuel ignites a propagating wave. To ensure that 
the details of the ignition process do not affect the later 
flame propagation (for instance, by inducing significant 
pressure waves, or by 'overdriving' the flame), we ran 
simple tests of using different temperatures to ignite the 
flame. Results for a Le = oo and Le = I Arrhcnius flame 
are shown in Fig. [f)| where we vary the ignition tem- 
perature by a factor of four. The two flames are shown 
because they have distinct ignition mechanisms — the 
Le = oo flame ignites due to thermal diffusion, eventu- 
ally heating a fuel layer to the point where significant 
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Fig. 6. — Position versus time, for a spherically-expanding 
Arrhenius flame with Le = co (left) and Le = 1 (right). Ignition 
for the two Lewis numbers proceeds differently. For each flame, 
the region x £ [0, 10] was hot ash at temperatures (bottom to top) 
T = 0.64,0.96,1.28,1.92,2.56, (only 0.64, 0.96, and 1.28 for the 
Le = 1 flame) with densities set to attain pressure equilibrium with 
the cold fuel at p = 1.0, T = 0.1 in the rest of the 100 cm domain. 
The adiabatic flame temperature is T;, = 0.64. As can be seen, 
the ignition times can vary significantly in the Le = co case (with 
the lowest-temperature case never igniting a propagating flame at 
all), but once ignition has occurred, the flame propagates with the 
same velocity. For the Le = 1 flame, even the ignition time is 
nearly constant. 



burning occurs, whereas the Le = 1 flame ignites due to 
combustible fuel diffusing into the already hot ash and 
igniting. We see that changing the ignition temperature 
affects the time-to-ignition for these flames, but it does 
not affect later flame propagation. 

4.6. Measuring curvature strain rate and Markstein 

number 

As discussed in |21 we are interested in the total di- 
mcnsionlcss strain, which for our two spherical flames 
are 

2// 



K&t = Kag + Ka 



■c 



If "gas 



2 J± 



2 dr f /dt 
l fZ 5o— (22) 



'7 



S? 



since drf/dt = w gas + S®. We have described measuring 
Tf above, and drf/dt can be calculated by differencing 
as long as care is taken to have sufficiently time-sampled 
data. Sf can be measured from our planar flame simula- 
tions after we have settled into a steady state. We have 
described methods to measure I / , but instead we choose 
to simply calculate Ka-r/^/, leaving 



S, = S?(1 



(1 



(i/Ma) (K&r/J,)) . 
£xKar/l/) 



(23) 
(24) 



so that a physical dimensional quantity, the Markstein 
length (Cm = i/Ma) can be measured directly by fit- 
ting, and then a dimensionless Markstein number can be 
calculated using any chosen thermal width of the flame. 

5. RESULTS 

The goal of this study is to investigate the effect of 
curvature on the speed and structure of several types of 
flames. As described above, the effects define the Mark- 
stein relation, so the majority of our results consists of 
quantifying this relation for each of the flames we con- 
sider. Best fits were calculated using a maximum likeli- 
hood linear least squares fit, and errors quoted in that fit 
come from variances of th e fit parameters unless stated 
otherwise (see for instance lBevingtonlll969|) . Astrophys- 
ical flame results are presented here; model flame results 
can be found in Appendix IO 

The planar flames were run first to determine the flame 
width and speed accurately. Experimentation showed 





Fig. 7. — Flame speed versus Kar/'/ for a 2.5 X 10 7 g cm~ 3 
(top left), 5 X 10 7 g cm -3 (top right), 1 X 10 8 g cm" 3 (bottom 
left), and 2 X 10 s g cm -3 (bottom right) pure carbon flame. 
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Fig. 8. — Flame speed versus Kay/i^ for a 5 X 10 7 g cm 
(left), 1 X 10 8 g i^" 1 3 ffem+m'} ■ i < 1 o v in^ it : ■>,, 3 
0.5 O flames. 



that we need at least 10 computational zones in the 10/90 
thermal width. We ran all of these flames with about 20 
points in this thermal width (as determined from the 
planar flame runs). Several densities were run, both in 
pure carbon and half carbon/half oxygen, and the flames 
were run both inward and outward in radius. Data on 
the planar flames is given in Tabled 

Fig. shows the flame velocity vs. scaled dimension- 
less strain for flame in a pure carbon medium, and Fig. [5] 
show results for a 50/50 Carbon/Oxygen medium. The 
results are summarized in Table Quoted errors on the 
Markstein length come from uncertainties in the fit. As 
is the case with the model flames, and is described in Ap- 
pendix El the data for the inward propagating flames is 
far more noisy than the outward propagating flames due 
to pressure waves from ignition transients being trapped 
inside the unburncd region. The noise is, however, rea- 
sonably symmetrically distributed about the line fit. We 
do not need the inward data to get a Markstein length, 
as the outward data alone is enough to provide this, 
so we can do fits of just the outward data to see how 
it compares to the complete data set. In the most se- 
vere cases, the difference in the Markstein length com- 
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Fig. 9. — Fits of the 1 X 10 s g cm -3 0.5/0.5 Carbon/Oxygen 
(left) and pure carbon (right) flames using only the outward prop- 
agating flame data. 



puted with and without the inward propagating data is 
10%. For the 1 x 10 8 g cm -3 pure carbon flame, Cm = 
— (8.366±.238) x 10~ 4 cm when all of the data is used, and 
Cm = — (9.003 ± .12) x 10~ 4 cm when only the outward 
propagating data is used. This flame had one of the the 
largest errors. The same density, 50/50 Carbon/Oxygen 
flame has C M = -(1.028 ± .003) x 10 -2 cm when all the 
data is used, and Cm = -(1.047±.002) x 10~ 2 cm when 
the outward propagating only data is used. Figure 
shows the fits for these two flame when the outward data 
only is used. They can be compared to Figures and 
In the chemical combustion literature cited here, mea- 
surements are often only made of outward- propaga ting 
flames (eg., Hassa n et al.lll998l IK won" fc Faethll2r)0ll) . 

To assess the effects of the boundary conditions on the 
results, we ran some flames with a domain twice the size 
as that used in the above results. Fig.llOlshows the speed 
vs. Ka/Z/ 1 x 10 8 g cm -3 pure carbon flame, where the 
domain size was 7.68 cm, compared to 3.84 cm used in 
the main study (Fig. 0|. The Markst cin length computed 
for this run is Cm = —(8. 159 ±.248) x 10 _4 cm, compared 
with the smaller domain value of Cm = — (8.366±.238) x 
10~ 4 cm. We see that these results are consistent within 
their fit errors. 

The tables show that the Markstein length increases 
as the density or carbon fraction decrease. The reason 
for this is the predicted dependence on Zeldovich num- 
ber. Table gives the Zeldovich numbers for each of the 
CF88 flames, found by fitting an Arrhenius rate to the 
measured e„uc from the simulation. A very good Arrhe- 
nius fit can be found for each flame, but the parameters 
vary from flame to flame. T a , for instance, varies slightly 
but only by about 10%. If we assume that T a is roughly 
constant and remember T5 3> T u for these flames (sec 
for instance Fig.0), then Ze « T a /Tb, so that increasing 
Tf, decreases Ze. 7& can be increased by increasing the 
fuel (Carbon) fraction so more energy is released or by 
increasing the density of the fuel and therefore, its de- 
generacy. (The more degenerate the fuel, the less of the 
energy released by burning will go into PdV work, and 
so the temperature of the ashes must increase instead). 
We see that the measured Zeldovich numbers do indeed 
follow these trends. 

The relationship between Ze and Ma is further ex- 
plored in Fig. ^2 We see that there is a trend of increas- 
ingly negative Markstein number with Zeldovich number, 
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Fig. 10. — Results of the 1 X 10 s g cm -3 pure carbon flame run 
with a domain twice the size used in the main study. This can be 
compared with Fig. 



Fig. 11. — Plot of the calculated Zeldovich number for the 
six astrophysical flames considered here and their measured Mark- 
stein numbers. The vertical error bars reflect the values Ma''' and 
Mat"', with the symbols marking the average of the two for a 
given flame. The bars with a star symbol ('*') represent the data 
from the 50% carbon flames, and those with a diamond symbol 
('o') represent the 100% carbon flames. Horizontal error bars rep- 
resent the uncertainty in measuring the Zeldovich number from 
fitting the burning rate. 



but there is at least one other variable involved as well; 
different carbon fractions give different trends. Also, as 
degeneracy lifts (for the p = 2.5 x 10 7 g cm -3 pure car- 
bon flame), the relationship grows more complicated. 

6. CONCLUSIONS 

We have reported on the curvature behavior of astro- 
physical and model flames. The behavior of astrophysical 
flames can be understood in terms of simple geometry; a 
flame bent outwards towards fuel will burn more slowly, 
because the diffusive heat transport is 'diluted', whereas 
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a flame bent inwards towards ash burns more quickly, 
because heat transport is concentrated. Quantitatively, 
we have seen that the astrophysical flames reported on 
have a reaction to a range of curvature- and flow-induced 
strain that is linearly related to the magnitude of the 
strain for the small strains (Ka^ < 1) examined here 
and have a Markstein length of size comparable to the 
flame thickness. 

Simplified model flames can have very different behav- 
iors than the astrophysical flames, if their structure is 
very different. In particular, a peaked burning rate and 
a separation between a preheat zone and a burning zone 
are responsible for astrophysical flame behavior on the 
scales explored here. To measure small-scale flame re- 
sponses, then, real fully resolved flames must be calcu- 
lated, or models which take into account strain effects 
ijMatalon fc Matkowskvl 11982ft must be used. Methods 
which use advect ion-diffusion equ ations but do not fully 
resolve the flame ljKhokhlovll995ft may work well at large 
scales, but would be inappropriate to use for modelling 
small-scale flame behavior. 

We have demonstrated that to measure accurate burn- 
ing velocities one must consider the bulk burning rate, 
rather than using differencing on the flame's position. 
It is difficult to accurately determine and subtract the 
fuel velocity, and relatively easy to accurately compute 
integral quantities. Even measuring the flame position 
required care, as small uncertainties in the measured po- 
sition can result in significantly different computed Mark- 
stein numbers. 

The magnitude of the response to strain increases with 
increasing Zcldovich number. For the flames examined 
here, decreasing carbon fraction or decreasing density 
both increase the Zcldovich number. However, the abso- 
lute magnitude of the response seen here has a stronger 
dependence on composition of the fuel than through the 
Zcldovich number. The magnitude of the response is 
also smaller than would be implied by the basic theory 
described in More complete derivations, including 
the decr ease in density a cross the flame, have been de- 
scribed ijSun et al.lll999l) . but the density jumps in the 
flames reported here are quite modest. Other fluid prop- 
erties (temperature, conductivity, diffusivity) change a 
great deal across the flame, more than the density jump 
would imply, because of the degenerate EOS; however, 
most theory in the terrestrial combustion literature as- 
sumes an ideal gas EOS. These results, then, must be ex- 
tended to the degenerate EOS applicable to white dwarf 
interiors. 

The Markstein lengths measured here represent the 
smallest scale wrinkling one can expect to see in astro- 
physical thermonuclear flames. Thus, it marks the lower 
scale end of the fractal-like behavior one might expect to 
see in large-s cale Landau-Darrieus instability growth, as 
described in ijBlinnikov fc Sasorovlll99fit iBlinnikov et al.1 
Il995l) . 

A velocity law such as in Eq. implies a burning ve- 
locity of zero — that is, quenching — for Ka = — 1/Ma 
or K = —S^/Cm- Clearly, extrapolating the behavior 
of the flames from the strain rates investigated in this 
paper to strains sufficient to quench the flame is prob- 
lematic, but keeping this in mind, we can use this as an 
estimate for the sorts of strain rates which can signifi- 
cantly affect local burning rates. A periodic sinusoidal 



shear flow of amplitude A and wavelength of n flame 
thicknesses has a peak strain rate of (2ir A) / '{nl /) . For 
that to be equal to the estimated quenching strain above, 
(A/S®) = n/(27rMa). Since the flames investigated had 
Ma « 1, we see that periodic strains on the flame thick- 
ness scales with velocities comparable to the flame speed 
can be expected to sig nificantly alter the flame's burning 
rate. See for example. iZingald l)2001ft . 

The direction of the astrophysical flame response to 
strain is to act against the strain, thus partially stabiliz- 
ing an initially planar flame against wrinkling by strain 
or instabilities. Given the measurements presented here, 
we can predict the flame's behavior under the Landau- 
Darrieus instability. From Eq. ^| we can calculate the 
change in the growth rate of the instability, and the 
largest mode stabilized, by the non-constant flame ve- 
locity. If, as in |21 k c = (a — l)/(2aZ/|Ma|) is the critical 
wavenumber for stability, we can calculate the largest 
wavelength stabilized, A c = 2ir/k c . These arc listed, for 
the astrophysical flames studied here, in Tablc|3J Results 
are given in terms of flame thicknesses (we use method 
I for concreteness) and in centimeters. We note that 
the increase in Markstein number as one goes to lower 
density approximately cancels the expected decrease in 
stability due to increasing density contrast as degeneracy 
lessens. 

Along with the growth rate of the Landau-Darrieus 
instability, the curvature behavior of the flame also 
affects its non-linear stabilization. As is well known 
l)Zeldovich et al.ll985f) . the growth of a perturbation in a 
propagating interface may be limited by cusp formation. 
This is not a property of the Landau-Darrieus instability, 
but simply of a propag ating interfac e by Huygens's prin- 
ciple (see, for example, IHechtlll98^) and can be demon- 
strated with the s o-called 'G-equation' (see for instance 
lLaw fc Sunel 120011) . There ha s been some astrop hysical 
interest in this very recently (|Ropke et alJl20f)3ti . Since 
cusp formation involves large curvatures, is very strongly 
affected by the curvature behavior of the flame, and in- 
deed, cusp formation can be completely inhibited by the 
flame behavior, as shown in Fig. ^| 

The local flame response to shear measured here makes 
it harder to significantly accelerate the burning on scales 
of a few tens of flame thicknesses by instabilities, and 
may not strongly affect larger-scale wrinkling. The same 
local flame response may, however, make a transition to 
detonation easier, as the burning rate changes signifi- 
cantly with shears of order /I®. 

This work has focused on single-step reaction mecha- 
nisms. Realistic flames, both astrophysical and terres- 
trial, involve many reactions and intermediate species. 
Future work will extend the results presented here to 
such multi-stage reaction networks. 
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Fig. A13. — On the left is the density profile of a spherically expanding Sedov explosion 0.05 seconds after the explosion at r = 0. The 
explosion energy was 1 erg, into an ambient pressure of 10 — 5 dynes cm -2 and an ambient density of 1 g cm -3 . On the right is the shock 
position of the same simulation over time. Simulated results are shown with open circles, and predicted results with a solid line. 




Fig. 12.— Solutions to the 'G equation', dG/dt = |VG|(1 + 
Cj^d 2 G/dx 2 ) tracking a propagating surface through a domain 
with no flow. The interface is originally a sinusoid of dimcnsionlcss 
amplitude 1, and wavelength 2n, and propagates downwards. If the 
flame propagates at a constant velocity of 1, it cusps as shown in 
the top left frame (at times 0, 1/3, 2/3, 1, 4/3). If the interface has 
a negative Markstein length equal to 1/20, the results are shown in 
the top right (at times 0, 1/3, 2/3, 1, 2, 3) and we see that cusping 
is postponed. Doubling the Markstein length gives behavior as 
shown on the bottom left, at times 0,1/3,2/3,1,2,3. Cusping is 
suppressed, and the flame begins to flatten. On the other hand, a 
positive Markstein length of 1/20 enhances cusping as shown on the 
bottom right (shown at times 0,1/3,2/3,1.) Numerical solutions 
of the PDE were calculated using Mathematica. 



simulations were performed with the Flash code, 
version 2.2. Some calculations were performed on the 
UCSC UpsAnd cluster supported by an NSF MRI grant 
AST-0079757. The astrophysical flame models pre- 
sented in the appendix arc available for comparison at 
http : //www.ucolick. org/~zingale/f lame_models/ 
The Flash code is freely available at 
http : //flash . uchicago . e du/| 



APPENDIX 
VERIFICATION TESTS 

In this paper, we describe results from the Flash code with additional modules for modelling diffusion and conduc- 
tivity, and for computin g hydrodynamics in 1-d sp herical coordinates. The Flash code has been well tested elsewhere 
without these modules ( Calder 2002; Frvxell 2000); in this appendix we present test results of the new modules added 
for the simulation of flames. 

1-D Spherical Hydro 

Since most of the simulations presented here were performed in spherical coordinates, we first test our solutions of the 
equations of hydrodynamics in spherical coordinates. The hydrodynamic algorithm and implementation used in the 
Flash code have been well-tested in Cartesian coordinates, so we need only test that the area factors and adaptive mesh 
refinement (AMR) routines are cor rectly implem ented for spherical coordinates in the Flash code. Shown in Fig. lA13l 
is the result of a Sedov explosion ljSedovlll959|) propagating in a 1-d spherically-symmetric domain. This calculation 
included AMR. The calculated result is shown with circles, and the predicted density profile ijLandau fc Lifshit3ll987j) 
is shown with a solid line. Also plotted is the shock position versus time. We see that the shock is well modeled with 
the Flash code in spherical coordinates. 



Flame Response to Curvature and Strain 



Unburncd state 






Unburned 


Burned 


P 




Ze 


a 


P (dyn cm 2 ) 


1 (IX) 


i itiyii cm ) 


1 (IS.) 


2.5 x 10' 


1 


13.4 


2.17 


3.14 x 10 24 


1.0 x 10' 


3.14 x 10 24 


4.05 x 10 a 


5 x 10' 


1/2 


13.5 


1.47 


8.27 x 10 24 
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Table 1 

Flame properties of the planar CF88 flames 



P ( g cm' 3 ) 


X c Sf (cm s" 1 ) 


l u ^> (cm) 


t\ [u) (cm) 


Cm (cm) 




Ma« 




2.5 x 10' 


1.0 


1.71 x 10 5 


6.53 x 10" 


-2 


3.86 x 10" 


-2 


-(3.107 ± .010) x 10" 


-'2 


-0.476 


-0.805 


5 x 10' 


0.5 


9.11 x 10 4 


5.28 x 10" 


- : 2 


3.07 x 10" 


-2 


-(5.496 ± .007) x 10" 


-2 


-1.04 


-1.79 


5 x 10 7 


1.0 


6.81 x 10 5 


9.01 x 10" 


-3 


5.49 x 10" 


-3 


-(4.837 ± .038) x 10" 


-3 


-0.537 


-0.881 


1 x 10 H 


0.5 


2.92 x 10 a 


9.93 x 10" 




6.84 x 10" 


-'i 


-(1.028 ± .003) x 10" 


-_' 


-1.03 


-1.50 


1 x 10 8 


1.0 


2.13 x 10 6 


1.85 x 10" 


-3 


1.31 x 10" 


-3 


-(8.366 ±.238) x 10" 


-4 


-0.452 


-0.639 


2 x 10 H 


0.5 


8.19 x 10 a 


2.33 x 10" 


-a 


1.72 x 10" 


-:i 


-(2.180 ± .014) x 10" 


-:•> 


-0.935 


-1.27 


2 x 10 s 


1.0 


5.56 x 10 6 


4.63 x 10" 


-4 


3.74 x 10" 


-4 


-(9.073 ± 2.44) x 10" 


-5 


-0.196 


-0.243 



Table 2 
cf88 flame data. 



p ( g cm- 3 ) X C \c/l° f {1) A c (cm) 



2.5 x 10' 


1.0 


11.1 (7.24 ±.02) xlO" 


l 


5 x 10' 


0.5 


50.8 (2.15 ±. 003) xlO u 




5 x 10 7 


1.0 


14.3 (1.28 ±.01 ) xlO" 


i 


1 x 10 s 


0.5 


48.9 (4.86 ±.01 ) xlO" 


l 


1 x 10 8 


1.0 


13.9 (2.57±.07 ) xlO" 


2 


2 x 10 M 


0.5 


52.3 (1.21 ±.007) X 10" 


1 


2 x 10 8 


1.0 


8.46 (3.91 ±1.07) x 10~ 


3 



Table 3 

Predicted critical wavelength for Landau- Darrieus 
instability, in terms of flame thickness (method i) 
and in centimeters, for the cf88 flames described 
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Fig. A14. — On the left, the evolution of the best-fit Gaussian widths as a function of time in the simulation, both predicted (solid line) 
and measured ('o'), for a constant diffusivity of D t ^ = 1 cm 2 s — 1 . The original profile is a Gaussian of width a = 1/32 cm; the resolution 
is Ax = cto/16 = 1/512 cm. Plotted on the right is the relative difference between the best-fit Gaussian width and the predicted Gaussian 
width. 



Diffusion Tests 

To test the diffusion module described in i)3.3l we diffuse a Gaussian pulse of temperature through a constant-density 
domain. As is well known, a planar initially Gaussian pulse will remain a Gaussian when diffused with a constant 
diffusivity D t ^, with a width 

(Al) 



<r(i) = y/a% + 2Dth(t-t ), 



where <7q = o~(t = to). We use this to test the thermal diffusion in Cartesian coordinates. An initial temperature pulse 
is created and evolved in the Flash code. For this test, hydrodynamics is not calculated, so that we can get a clean 
test of the diffusion. 

Quantitative results arc shown in Fig. I AMI The measured Gaussian widths were computed by fitting the evolved 
profile with a Gaussian of arbitrary amplitude and width. As we see, we get excellent agreement over a diffusive 
timescale. The same test was run with a species diffusion problem; identical results were obtained, since the code for 
the two diffusion operators is identical. 

A 1-d spherically symmetric Gaussian profile is not a solution of the spherically-symmetric diffusion equation 

dT 
~dt 

Instead, the solutions are spherical Bessel functions. Thus we consider a temperature profile 



I d 2 d 



(A2) 



T(r,t)=T (l+j (kr)e 



-k A D ttl t 



) 



(A3) 



We simulated a diffusing Bessel function profile with k = 1 cm -1 , To = 1. K, and D t h = 100 cm 2 s — 1 . A clo se-up 
of the calculated and predicted profiles and the error after 2 diffusion times (0.02 seconds) are shown in Fig. IA15I 
The error is small everywhere except at the right boundary, because the infinite domain in which the spherical Bessel 
functions form a solution to the diffusion equation is modeled by a truncated domain with a zero-gradient boundary 
condition. Reflecting boundary conditions are used at the left boundary. 

Boundary Conditions 

Verifying a boundary condition is difficult; here we simply present evidence that the boundary condition allows the 
material from the flame to smoothly leave the domain. Fig. IA161 shows a space-time diagram of points in the domain 
of a planar flame simulation described in 21 The Lagrangian trajectories were calculated as a post-processing step on 
the results of the simulations, taking a set of initial points and tracking them through the velocity field. When a point 
is overtaken by the flame, it expands and falls behind the flame, eventually leaving the domain; material ahead of the 
flame moves slowly in the other direction. In both cases, the material must advect smoothly out of the domain. We see 
in Fig. lAlfil that material advects out of the domain without any artifacts. The initial conditions for this simulation 
was hot (T = 0.96) ash on the left up to x = 10, cold (T = 0.1, p = 1) fuel on the right, in pressure equilibrium, and 
Lc = oo. A Zc ~ 5 Arrhcnius burning rate and a 7 = 5/3 ideal gas EOS were used. Ignition occurred after about 70 
seconds. 



Flame Speeds and Structure 
Model Flames 

For our model flames, in Table lA"4l we summarize the results of a resolution study for the Le = 00, Ze « 5 flame. We 
see consistent resolution requirements for both these flames and the astrophysical flames listed above. We use a finest 
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Fig. A15. — On the left, the temperature perturbation after two diffusion times ('o'), and the predicted profile (solid line). On the right, 
the error in the evolved temperature perturbation at the same time. Error is small everywhere except near the boundary conditions. 




position (cm) 



Fig. Af6. — Points moving along Lagrangian trajectories as a planar Arrhenius model flame ignites at x = fO at time t fa 70. Material 
on both sides of the flame smoothly leaves the domain. 



resolution of 1.56 x 10 -2 cm, which is seen to provide well-converged results, for the Zc ~ 5 Arrhenius flames and the 
KPP flames described here, even though the KPP flames and finite Lewis-number Arrhenius flames are thicker than 
the infinite Lewis number Arrhenius flame. 

CF88 flames 

The properties of astrophysical carbon flames were investigated in great detail in lTimmes fc Wooslevl (j!992fl . There, 
several different methods and nuclear reaction networks were tried. All of the reaction networks they considered in- 
volved multiple reactions. The absence of multiple reactions in our study makes the temperature structure considerably 
simpler. After the initial carbon burns, the temperature remains flat. With a more extensive network, later reactions 
would continue to heat (or possibly cool) the region behind the flame. The 10/90 flame width definition would be 
more sensitive to these later reactions, and would tend to produce a wider flame than those we get. 

Flames were simulated propagating through a range of densities relevant to white dwarf interiors. Of particular 
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^zoncs 

1.25 x 1CT 1 7 4.046 X 10~ 2 

6.25 x 10~ 2 14 4.001 x 10~ 2 

3.13 x 10~ 2 28 3.999 x 10~ 2 

1.56 x 10- 2 57 3.999 x 10~ 2 



Table A4 

Numerically determined laminar planar flame velocity as a function of resolution in the Flash code. The 
resolution is expressed as the zone size, and the number of zones across the thermal thickness (iy^') of the 

Le = oo, Ze ~ 5 Arrhenius flame. 



Sx (cm) n Z oncs cms 1 /^"^ l ^ 11 ' 

5.21 x 10~ b 34 2.135 x 10 b 1.75 x 1CT 2 1.22 x 10~ 2 

1.04 x 10" 4 18 2.132 x 10 5 1.85 x 10~ 2 1.31 x KT 2 

2.08 x 10" 4 10 2.162 x 10 5 2.04 x 10" 2 1.51 x TO" 2 

4.17 x 10~ 4 6 2.194 x 10 5 2.47 x 10" 2 1.64 x 10~ 2 



Table A 5 

Numerically determined laminar planar flame velocity for the 1 x 10 8 g cm -3 pure carbon CF88 flame as a 
function of resolution. The resolution is expressed as the zone size and the number of zones across the 
thermal thickness of the flame. We see that the speed converges at ROUGHLY 20 POINTS in the laminar flame 

thickness. 



interest is determining how much spatial resolution is required to accurately simulate a flame. The results presented 
here were u sed as the benchmark for the flame propagation experiments presented in this paper. 

Fig. IA171 shows the structure of a 5 x 10 7 g cm pure carbon flame, propagating to the left. We see that the large 
jump in temperature happens right at the peak nuclear energy generation rate, and then the temperature is almost 
perfectly flat behind that. The temperature jump is accompanied by almost a factor of two density drop behind the 
flame. 

An accurate description of the deflagration speed requires resolving the front. For all the flames studied here, we 
ran with several different resolutions, in order to find out how much resolution was necessary for the flame speed and 
width to converge. Table lA~5l summarizes the results of one such resolution study, for the 1 x 10 8 g cm -3 pure carbon 
flame. We see that 20-40 zones in the flame width is sufficient to determine the flame speed accurately. The main 
calculations in this paper were run with a resolution that put about 20 computational zones in the flame width (1°^). 
We did not consider simulations with less than 5 zones in the flame width. 



conservative interpolation 

In FLASH, the computational zones are organized into blocks, typically containing 8-16 uniformly gridded zones 
in each coordinate direction. Adapting to the flow is achieved by allowing blocks to vary in resolution by a factor of 
two with respect to neighboring blocks. When a region needs to be refined, the parent block is cut in half in each 
dimension, creating 2 d new blocks, where d is the dimensionality. Initializing the newly created zones in each child 
block requires that we reconstruct the zone-average data in the parent (and its neighbors) and use this to compute the 
child data in a conservative fashion. This method is also needed at jumps in resolution, when exchanging guardcell 
information between neighboring blocks. 

Reconstruction is achieved by fitting a quadratic polynomial to the parent data. Since Flash a finite-volume code, 
the coefficients of the in terpolation polynomial must be chosen to reproduce the zone average values of the parents in 
the stencil. Figure IB 1 81 shows the three coarse zones used to initialize the children. 

Given the polynomial, 

f(r) =a ir 2 +a 2 r + a 3 , (Bl) 

we apply the constraints, 

Air f ri + 1 / 2 

tt / rHrf{r) = (f) l , (B2) 

to each of the parent zones and solve to find a±, a 2 , and 03. Here, Vi is the volume of the zone, which, for spherical 
coordinates, is 

Vi = ^(rf +1/2 - r 3 _ 1/2 ) , (B3) 

and (f) i is the average value of variable / in zone i. The zone edges arc denoted by half-integer indices, as illustrated 
in Figure IB 181 This polynomial can then be integrated over the children to initialize their data. It is also used when 
filling guardcells at coarse-fine interfaces. 
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FlG. A17. — Flame structure for a planar 5 x 10 7 g cm 3 pure carbon CF88 flame. Shown are temperature (top left), density (top 
right), carbon mass fraction (bottom left), and nuclear energy generation rate (bottom right). 
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Fig. B18. — Schematic showing three coarse zones with zone-averaged data (f) i _ 1 , (f)n and (f) i _^ 1 (top). The i-th zone is refined to 
create two children (shown on bottom), whose data is to be initialized with (f) ic and (/} ic+1 , through conservative interpolation. 



Applying Eq. IB2l to the i-th zone yields 



^aiQln + ^0803^ + as = {f)i (B4) 



where, for notational convenience, we define the operator 



e:^ r /2 - r 1/2 m 



r i+l/2 T i-l/2 
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We also use the i — 1 and i + 1 zones, and get 

>5 



^aiGs + ^203 +a 3 = (/},;_! 
^ai03 r i+i + ^a2G3 r i+i + a 3 = {f) i+ i 



Equations IB4HB7I can be solved for ai, et2, and 03 to yield: 



a 2 





03 n ~ 03 r i+l 




03 r » - 03 r *-l 


{{f)i 


03 r * - 03 r »-l 


03 r i - 03 r i+l 




03 r » - 03 n+i 


03 r i - 03 r i-l 



dl 



(B6) 
(B7) 

(B8) 
(B9) 

g -™-. 4 -^3M- (BIO) 

Together with Eq. IB1I these reconstruct the data in the i-th zone. When we re fine, two children are created from 
the i-th zone, each taking up half the interval. This is illustrated in Figure ITj 181 on the lower axis. For convenience, 
we use ic to ref er t o the child index space, with r ic _!/ 2 = ^i-i/2) as shown in the figure. The child data is found by 
integrating Eq. IB1I over the child zone as 

47T r' c+1/2 o 

r 2 drf(r) , (Bll) 

ic-1/2 



03^ - 03 r i+l 



(f)i - - GS r * - 05 r «+i 



and 



fl 3 = (f)i - |ai03 r i - 7 fl 20^ r 



</) ic = 



V it 



where ic refers to the child zone index. Carrying out the integration, we find the value of the child data, in terms of 
the parent data and it's immediate neighbors is 



(f) tc = j|ai03 r ic + ^a 2 03 r lc + a 3 . 



(B12) 



This procedure is guaranteed to be conservative, but not guaranteed to be monotonic. If monotonicity is a require- 
ment, then we can check whether the newly initialized child data falls outside the limits of the coarse zones in the 
stencil, and if so, fall back to direct insertion. 

MODEL FLAMES 

In this appendix, we describe methods and results for 'model' flames, both to connect with the combustion literature 
and to provide simpler models for understanding the astrophysical flames. 

Equation Of State 

For the simple model flames we use a ideal gas gamma-law equation of state 

P = (7 - l)pe , 

supplemented with an expression relating the temperature and the pressure (ideal gas approximation), 



N a k 

= ~i~p T ' 



(Cl) 
(C2) 



where e is the internal energy / gram, iV a is Avogadro's number, k is Boltzmann's constant, and A is the average 
atomic mass of the mixture of fluids. For the model flames, we use 7 = 5/3. 

Diffusion 

For some of our model flames, we allow finite Le, and so species diffusion is also calculated. Here, the diffusion term 
of the species evolution equation is 

d P X l 



dt 



(C3) 



where Xi is the mass fraction of species i, and D s is the species diffusivity. When species diffusion is included, D s 
is assumed constant and fixed. The species diffusion is solved in the same way as the thermal diffusion, described in 

Burning 

In § 13.5.11 we described the burning rates for the astrophysical flames. Here we describe simpler burning rates used 
for the model flames. 
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Lc 


S° (cm s- 1 ) 


(cm) 


t\ [U) (cm) 


C M M (cm) Cm^ (cm) 


1 


0.0896 


5.43 


4.42 


-.533 ±.139 +1.28 ±.60 


oo 


0.0927 


5.31 


4.23 


-.604±.139 ±1.92±.82 



Table C6 

KPP FLAME DATA. 



KPP 

A simple burning mech anism, often used by applied mathematicians because it is tractable analytically, is due to 
iKolmogorov et al"1l)1937f) . known for the authors of that paper as KPP. It is a one-step reaction that irreversibly 
converts fuel to ash: 

d ®=- d 4l = aeX f (04) 
dt dt 1 y ' 

where is a linearly scaled temperature such that = 1 corresponds to the temperature of the burned ash, and 
9 = corresponds to the temperature of the unburned fuel. The concentration (mass fraction) of fuel is given by Xr , 
also varying from 0—1. An arbitrary parameter a sets the rate of reaction. 

Along with a and the diffusivities D t h and D s , we have another parameter we can choose — the energy input of the 
burning with respect to the background energy. We parameterize this by the ratio of the energy released by burning 
a gram of material to the thermal energy of a gram of fuel in the unburned state 

(C5) 

where Ae is the binding energy released by burning a single particle of fuel, \if its mass, and Cp is the heat capacity 
at constant pressure of the fuel. 

A rrhenius 

Real flames, both astrophysical and terrestrial, have such strongly peaked temperature dependencies that their 
structure can be approximately be broken up into a diffusive "preheat" zone, where reactions are negligible, and a 
thin reaction zone, where most of the burning occurs and the temperature is roughly constant. Model flames using a 
KPP burning rate lack this important feature, because of KPP's fairly gentle tem perature depend ence. 

Thus, as a more realistic model, we consider an Arrhenius rate (see for example lWilliamslll985|) 

^ = aT a p n X?e- T «/ T , (06) 

where T a is an 'activation temperature' large compared to any of the temperatures in the system. If the molecular 
weights of the fuel and ash are the same, then in terms of Xf this becomes 

^ = - a ^ rP n - 1 XTe- T "/ T . (07) 
dt qT u 

Results 

For the scale-free model flames, dimensional values of parameters are somewhat arbitrary. Parameters were chosen 
so that the fastest outward-propagating spherical flames would have drj /dt « 1 in a domain where the slowest sound 
speed would be c s = 100, and all flames were run with these parameters. D t h was fixed at 0.0625 cm 2 s _1 , and 
Dg = Dth/Le. The unburned state was set to T u = 0.1, p u = 1. and the energy release was set so that TjP = 0.64. The 
atomic weight of both the fuel and ash were set to 1385.1. 

We consider tw o sets of flames — 'KPP names', which have an ideal gas EOS and burn with the KPP burning 
rate describe d in 3C.3.1I and 'Arrhenius flames', which also use the ideal gas EOS but burn with an Arrhenius law, 
described in 30.3.21 Simulations were run with varying temperature dependence and Lewis numbers. Simulations with 
finite large Lewis numbers become computationally increasingly expensive, as an increasingly smaller species diffusion 
zone must be resolved adequately. 

KPP flames 

Results for KPP flames are shown in Fig lC19l Shown are results from simulations with a = 0.3, Le = 1, and a = 0.3, 
Le = oo. Data from outwardly-propagating and inwardly-propagating spherical flames, and planar flames are shown. 
Quantitative results are summarized in Table IC6l The resolution used was that chosen in the Appendix for the Ze ~ 5, 
Le = oo flame, although the KPP flames are much thicker. 

We use the e nuc -based method for flame velocities as described in 34.41 even though the approximation of a thin 
reaction front is a poor one for these KPP flames. The quantitative results shown here can be changed significa ntly 
by choosing the flame position differently within the burning zone of the flame; the errors quoted for Cm in Table IU6l 
come from using the approximate ±1 cm uncertainty in flame position shown in Fig. EI 
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Fig. C19. — Flame speed versus Kay/ij- for KPP flames, with Le = 1 (left) and Le = oo (right). For each graph, data points for positive 
KaT/lf are from simulations of a spherically expanding flame, data points for Ka^/Z/ = are for a planar flame, and those for negative 
Kay/ij- are for a spherical flame propagating inwards. 
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1.326 


0.5422 


-0.297 ±0.001 


-0.224 


-0.548 




5 


0.0326 




1.025 


0.4052 


-0.654 ±0.001 


-0.638 


-1.614 



Table C7 
Arrhenius flame data. 



We see here that KPP flames do not behave in a way described by the Markstcin relation. Although the negatively 
and positively strained flames separately respond significantly and linearly to strain, the sign of the response is different 
in the two cases. Straining the KPP flames here result in a slowing down of the flame, regardless of the sign of the 
strain. 

That the KPP flame responds differently to strain rates of order a flame crossing time than other flame models is 
easily understood. Where the astrophysical flames described in JHhave a distinct burning and preheat region, the KPP 
flames do not. Thus, in this case, the straining stretches not just the preheat zone — essentially the preconditioning 
for the burning zone — but the burning zone itself. 

Arrhenius flames 

Arrhenius flames were run with Ze ~ 5 and 8. The flame structures for Ze ~ 5 planar flames with different Lewis 
numbers are shown in Fig. IC20I Le = oo flames with positive, zero, and negative strain arc shown in Fig. IC21I The 
plots of flame speed versus scaled dimensionless strain are shown in Fig. IC 3221 and Fig. IC23I Table IC7I gives the 
quantitative results from these flames. Quoted errors on Cm come from uncertainties in the best fit. The Le = oo, 
Ze ~ 8 flame proved too difficult to reliably ignite with the constant diffusivity used here. 

The Arrhenius flames response to strain is linear, and the Markstein lengths are of order the flame thickness or 
smaller. The sign of the response is negative for Le = oo flames, as expected. The sign changes for smaller Lewis 
number as the effect of fuel diffusion into the ash, which acts in the opposite direction of the thermal diffusion into 
the fuel, becomes more significant. 

One sees in these results and in the astrophysical results that flame response to negative curvature is noisier than 
that to positive curvature. This is caused by the physical setup of the simulations. For the ingoing flames, pressure 
waves caused by ignition transients can be trapped inside the spherical flame due to the density jump at the flame's 
position. In our f-d simulations, the pressure waves bounce along the interior of the flame. The waves bounce off of 
the reflecting boundary condition, correctly modelling the same wave coming from the part of the spherical domain on 
the other s ide o f the origin. The pressure waves are only amplified as the flame move s inwa rds. This is shown on the 
left of Fig. IC24l for the Le = 5, Ze ~ 5 ingoing Arrhenius flame. On the right of Fig. IC24l is plotted the same values 
for the corresponding outward-propagating flame. Here, the transient pressure waves can largely leave the domain. 

The pressure fluctuations, while small (« 0.2% for the ingoing flames, and w 0.02% for the outgoing flames), are 
enough to slightly modify the local burning, leading to the observed scatter in measured burning velocities. As the 
temperature sensitivity of the flame increases, the noise caused by the same pressure fluctuations become larger; this 
explains the increased noise in the Ze ~ 8 flames. 
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Fig. C20. — Shown is the flame structure (solid line is scaled e n uc, dashed line is Xf, and dotted temperature) for planar Arrhcnius 
Ze ~ 5 flames, with Le = 1,2, 5, oo at the top left, top right, bottom left, and bottom right. e nuc is scaled to 81434.6 ergs g _1 s , and 
temperature to the adiabatic temperature of the Le = oo flame, T® = 0.64. As the species diffusion decreases, the size of the location where 
there is both enough heat and fuel to burn decreases, and the burning zone narrows considerably, increasing peak e nuc (although the total 
burning rate only changes by a factor of 2.) X-axis is distance into the unburned material from the flame position rr. 
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Fig. C21. — Shown is the flame structure (solid line is scaled e n uc, dashed line is Xf, and dotted temperature) for Arrhenius Ze ~ 5 
flames which are the outward propagating (top), planar (middle), and inward propagating (bottom). Ka^/i/ is +.2625 cm" 1 for the 
outward propagating flame, cm -1 for the planar flame, and —.03675 cm -1 for the inward flame. Although the change in the thermal 
and fuel structure is very small, the flame speed changes by ~ 10%. 
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Fig. C22. — Flame speed versus Kax/lf for Ze ~ 5 Arrhcnius flames. Shown arc plots for Le = 1 (top left), Le = 2 (top right), Le = 5 
(bottom left), and Le = oo (bottom right). For each graph, data points for positive Ka^/lf are from simulations of a spherically expanding 
flame, data points for Ka^/lf = are for a planar flame, and those for negative Ka^/lf are for a spherical flame propagating inwards. 
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Fig. C23. — Flame speed versus Ka^/Z/ for Ze ~ 8 Arrhenius flames. Shown are plots for Le = 1 (top left), Le = 2 (top right), and 
Le = 5 (bottom). For each graph, data points for positive Ka^/// are from simulations of a spherically expanding flame, data points for 
Kar/ij- = are for a planar flame, and those for negative Kay/// are for a spherical flame propagating inwards. 
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Fig. C24. — Scaled pressure fluctuation ((P — Pq)/Po) m the computational domain, with scaled e nuc to show the flame's position, for an 
inward-propagating Ze ~ 5, Le = 5 Arrhenius flame on the left, and the corresponding outward-propagating flame on the right. Pressure 
waves can be trapped inside the flame (eg., in the unburncd fuel) for the inward flame, and amplified as the flame moves inwards; these 
fluctuations in the fuel state lead to the observed scatter in burning rates. The pressure waves are largely able to leave the domain in the 
outward-propagating case, causing reduced scatter. 
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